!************************************************************
! block data of common parameters for the stratification 
! equations
!************************************************************
module commondat
 implicit none

 integer (kind=4), public :: l 
 integer (kind=4), public :: imr
 real (kind=8),    public :: x0,xmax
 real (kind=8),    public :: tt00, rho00
 real (kind=8),    public :: rg, cp, g
 real (kind=8),    public :: m1, m2, m3    !three layer index param. 
 real (kind=8),    public :: m0            !one layer index param.
 real (kind=8),    public :: d, xtac, xtop 
 real (kind=8),    public :: x1, x2, x3    ! four layer poly profile

 !GRID SETTINGS________________________
 !_____________________________________    (S01 - runs)
 data l             / 240 /
 data x0,xmax       / -0.2, 1.2 / 
 
 !INITIAL THERMODYNAMICS_______________
 !_____________________________________   !(S01 - runs)
 !data tt00, rho00   / 200., 500./
 !data rg, cp, g     /  0.4,   1., 140./ 
 data tt00, rho00   / 100., 200./         !(testing ...)
 data rg, cp, g     /  0.4,   1., 65./
 !data tt00, rho00   / 100., 200./        !(S2.01 - runs) 
 !data rg, cp, g     /  0.4,   1., 70./
 
 !POLYTROPIC INDEX PROFILE PARAMETERS__
 !_____________________________________    (S01 - runs)
 data m1, m2, m3    /   2., 1.49,  1.4/
 data m0            / 1.5 / 
 data d, xtac, xtop / 0.05,  0.0,  0.9/
 data x1, x2, x3    /  0.0, 0.85,  1.1/
 
 !FLAGS________________________________
 !_____________________________________
 !
 !   imr = n choose an n-layer poly-
 !           tropic profile.          
 !_____________________________________
 data imr           /  3  /  

 
end module commondat
!************************************************************
!************************************************************
! MAIN ...
!************************************************************
!************************************************************

program main

 use commondat
 implicit none

 call star ()

 stop
end

!************************************************************
!************************************************************

subroutine star ()
 
 use commondat
 implicit none

 integer (kind=4),parameter :: n=3
 integer (kind=4), parameter :: out_unit=20
 integer (kind=8), parameter :: out_unit2=10
 integer (kind=8), parameter :: out_unit3=30
 integer (kind=4) k
 external star_f
 real (kind=8) x, xaux, dx
 real (kind=8) rhiso0       !isoentropic base profile
 real (kind=8) u0(n),u1(n)
 real (kind=8) the0, the1
 real (kind=8) cap
 real (kind=8) ss
 real (kind=8) hs(n)
 real (kind=8) mr_out
 logical :: file_exist
 character (len=1024) :: cmd

 dx=(xmax-x0)/(l-1)
 cap=rg/cp
 u0(1)=tt00
 u0(2)=rho00
 u0(3)=rg*rho00*tt00
 the0=u0(1)
 ss=g/(cp*the0)

 INQUIRE(FILE="polyprof.txt", EXIST=file_exist)
 if(file_exist) then
  WRITE(*,*)  "files of a previous running already exists, deleting ... "
  cmd = 'rm data1.txt scales.txt polyprof.txt'
  call system(cmd)
 endif 
  
 WRITE(*,*)  "printing polytropic profile ..."
 
 !prints poly profile
 x=x0
 do k=1,l
  call polyprof(x,mr_out)
  open(unit=out_unit3,file="polyprof.txt",action="write", position="append")
  write(out_unit3,'(2x,g14.6,2x,g14.6)') x, mr_out
  close(out_unit3)
  xaux=x+dx
  x=xaux
 end do

 WRITE(*,*)  "calculating stratification profiles ..."
 
 !computes and prints the stratification profiles
 x=x0
 do k=1,l
  rhiso0=rho00*(1.-ss*(k-1)*dx)**((1./cap)-1) !isoentropic base profile
  open(unit=out_unit,file="data1.txt", action="write",position="append")
  write(out_unit,'(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)') x,u0(1),u0(2),the0,rhiso0
  close(out_unit)
 
  call scaleh(x,n,u0,hs)
  open(unit=out_unit2,file="scales.txt", action="write",position="append")
  write(out_unit2,'(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)') x, hs(1), hs(2), hs(3)
  close(out_unit2)
  
  xaux=x+dx
  call rk4vec(x,n,u0,dx,star_f,u1)
  the1=u1(1)*((tt00*rho00)/(u1(1)*u1(2)))**cap
  
  x=xaux
  u0(1:n)=u1(1:n)
  the0=the1
 end do

 WRITE(*,*)  "DONE ..."

 return
end

!************************************************************

subroutine star_f (x,n,u,uprime)

 use commondat
 implicit none

 integer (kind=4) n
 real (kind=8) x
 real (kind=8) u(n)
 real (kind=8) devs(n)
 real (kind=8) uprime(n)
 real (kind=8) mr
 
 call fielddevs(x,n,u,devs)
 uprime(1)=devs(1)
 uprime(2)=devs(2)
 uprime(3)=devs(3)

 return
end

!************************************************************

subroutine polyprof(x,mr)

 use commondat
 implicit none

 real (kind=8) x, mr
 
 if(imr==3) then
  mr=m1-(m1-m2)*0.5*(1+erf((x-xtac)/(d)))-(m2-m3)*0.5*(1+erf((x-xtop)/(d)))
 else if(imr==4) then
  mr=m1-(m1-m2)*0.5*(1+erf((x-x1)/(d)))-(m2-m3)*0.5*(1+erf((x-x2)/(d)))-(m3-m1)*0.5*(1+erf((x-x3)/(d)))
 else if(imr==1) then
  mr=m0
 endif

 return
end

!************************************************************

subroutine fielddevs (x,n,u,fdev)
  
 use commondat
 implicit none
 
 integer (kind=4) n
 real (kind=8) x
 real (kind=8) u(n)
 real (kind=8) fdev(3)
 real (kind=8) T, rho, P
 real (kind=8) mr

 call polyprof(x,mr)

 fdev(1)=-g/((mr+1)*rg)                 !dT/dr 
 fdev(2)=(u(2)/u(1))*(-g/rg-fdev(1))    !d(rho)/dr
 fdev(3)=-u(2)*g                        !dP/dr

 return
end

!************************************************************
! Computes the Temperature, density and pressure scale
! heights of the corresponding stratification profiles
!************************************************************
subroutine scaleh (x,n,u,scales)
 
 use commondat
 implicit none 

 integer (kind=4) i, n
 real (kind=8) x
 real (kind=8) u(n), devs(n), scales(n)

 call fielddevs(x,n,u,devs)
 do i=1,n
  scales(i)=-1./(devs(i)/u(i))
 end do
 
 return
end 

!************************************************************
! Fourth order Runge-Kutta 
!************************************************************
subroutine rk4vec (t0,m,u0,dt,f,u)

 implicit none

 integer (kind=4) m
 real(kind=8) dt
 external f
 real (kind=8) f0(m),f1(m),f2(m),f3(m)
 real (kind=8) t0,t1,t2,t3
 real (kind=8) u(m),u0(m),u1(m),u2(m),u3(m)

 call f (t0,m,u0,f0)

 t1=t0+dt/2.0D+00
 u1(1:m)=u0(1:m)+dt*f0(1:m)/2.0D+00
 call f (t1,m,u1,f1)
 t2=t0+dt/2.0D+00
 u2(1:m)=u0(1:m)+dt*f1(1:m)/2.0D+00
 call f (t2,m,u2,f2)

 t3=t0+dt
 u3(1:m)=u0(1:m)+dt*f2(1:m)
 call f (t3,m,u3,f3)

 u(1:m)=u0(1:m)+dt*(f0(1:m)+2.0D+00*f1(1:m)+2.0D+00*f2(1:m)+f3(1:m))/6.0D+00
!debugging the rk4 coefficients
! print*, t0,u1(1:m),u2(1:m),u3(1:m) 

 return
end
